% DO THIS BEFORE EXECUTING
addpath(genpath('/media/sda2/Users/Nermine/Desktop/code/sdp/gurobi301/gurobi_mex_v1.20'));

clc;
clear all;
close all;

shouldwePlot = 1;
% Generating Data
Ntrain = 1800; % training will contain Ntrain examples.
l = [ones(Ntrain/2,1); -ones(Ntrain/2,1)];              % label #Ntrain
d = [l/2 + randn(Ntrain,1)/1.5  l/2 + randn(Ntrain,1)/1.5]; % data #Ntrain

triangledataflag=1;
if(triangledataflag==1)
    %or generate two traingular data clusters to illustrate the point
    %get a square
    xrandgen1 = rand(Ntrain,1);
    xrandgen2 = rand(Ntrain,1);
    xrandgen = [xrandgen1 xrandgen2 ones(Ntrain,1)];
    poslowerTri = find(xrandgen(:,1)<xrandgen(:,2));
    posupperTri = find(xrandgen(:,1)>=xrandgen(:,2));
    %translate the points appropriately
    T = [1 0 0.99;
         0 1 -0.99;
         0 0 1];
    xpositionedpart = xrandgen(poslowerTri,:)*T';
    % plot(xpositionedpart(:,1),xpositionedpart(:,2),'b.'); hold on;
    % plot(xrandgen(posupperTri,1),xrandgen(posupperTri,2),'g.'); hold off;
    xpositionfull = [xpositionedpart;xrandgen(posupperTri,:)];
    xpositionfull = 3*xpositionfull*[cos(90*pi/180) -sin(90*pi/180) 0; sin(90*pi/180) cos(90*pi/180) 0; 0 0 1 ]';
    d = xpositionfull*[1 0 0; 0 1 -1; 0 0 1 ]';
    d = d(:,1:2);
    %plot(d(:,1),d(:,2),'.');
    l  = [ones(length(xpositionedpart),1); -ones(Ntrain-length(xpositionedpart),1)];
    flipprob = rand(Ntrain,1);
    for i=1:length(flipprob)
        if(flipprob(i)>0.9) 
            l(i) = l(i)*-1;
        end
    end
end

trainingdata = [d(1:Ntrain,:) l(1:Ntrain)];             % appended #Ntrain

B = glmfit(trainingdata(:,1:2), [0.5*trainingdata(:,3)+0.5 ones(length(trainingdata(:,3)),1)], 'binomial', 'link', 'logit');
Ftrain=trainingdata(:,1:2)*B(2:end) + B(1);   %%%% This is the model
%legacy start
Proba1=ones(length(trainingdata(:,3)),1)./(ones(length(trainingdata(:,3)),1)+exp(-Ftrain));
Proba0=ones(length(trainingdata(:,3)),1)-Proba1;
% Finding pts whose predicted labels are 1,0 & between threshold centered at 0.5 resp
thresh=0.65;
Index0training=find(Proba0>thresh);
Index1training=find(Proba1>thresh);
IndexErrorZonetraining=find(abs(max(Proba1,Proba0)-0.5)<thresh-0.5);
% Plotting training data only with the decision boundary.
if(shouldwePlot == 1)
    figure;
    scatter(trainingdata(Index0training,1),trainingdata(Index0training,2),'r.'); hold on;
    scatter(trainingdata(Index1training,1),trainingdata(Index1training,2),'b.');
    scatter(trainingdata(IndexErrorZonetraining,1),trainingdata(IndexErrorZonetraining,2),'g');   hold off;
    axis equal;
end
%legacy end
trainingloss = sum(log(1+exp(-(trainingdata(:,3).*Ftrain))));

%first 5 test points
lt = [1; 1; 1; -1; -1];
dt = [lt/2 + randn(5,1)/5  lt/2 + randn(5,1)/5];
testdata = [dt lt];
%the 6th point
xt = [-1:0.1:1];                          epsilonval = 1e-6;
m = -(B(2)/(B(3) + epsilonval));    c = -B(1)/(B(3)+epsilonval);                     
yt = m*xt +c;
testdata(6,:) = [xt(end) yt(end) -1];

% Plot of the training and the 5 points.
if(shouldwePlot==1)
    figure;
    pos = find(trainingdata(:,3)==1);
    plot(trainingdata(pos,1),trainingdata(pos,2),'b.');
    pos = find(trainingdata(:,3)==-1);  hold on; 
    plot(trainingdata(pos, 1),trainingdata(pos, 2),'r.');
    pos = find(testdata(:,3)==1);
    plot(testdata(pos,1),testdata(pos,2),'yo');
    pos = find(testdata(:,3)==-1);  hold on; 
    plot(testdata(pos, 1),testdata(pos, 2),'mo');
    plot(xt,yt);
    axis([-3 3 -3 3]);                  hold off;
end


%Find probabilities q on test data (of size 6)
Ftest=testdata(:,1:2)*B(2:end) + B(1);   %%%% This is the model
q=ones(length(testdata(:,3)),1)./(ones(length(testdata(:,3)),1)+exp(-Ftest));

%Plot to visualize the probability surface. needs xt vector from above.
% x2 = -1:0.1:1;
% for i=1:length(xt)
%     for j=1:length(x2)
%         psurf(i,j) = 1./(1+exp(-(xt(i)*B(2) + x2(j)*B(end) + B(1))));
%     end
% end
% figure; surf(psurf)

%for legacy.
pos_test = [1 2 3 4 5 6];
pos_l1 = [1 2 3];
pos_l0 = [4 5];
pos_t6 = 6;

%legacy start
%plotting the test nodes
if(shouldwePlot == 1)
    figure;
    scatter(testdata(pos_l1,1),testdata(pos_l1,2),'r'); hold on;
    scatter(testdata(pos_l0,1),testdata(pos_l0,2),'b');
    scatter(testdata(pos_t6,1),testdata(pos_t6,2),'c^'); axis([-3 3 -3 3]);
    hold off;
end

%checking probabilties of test nodes
%Assigning probabilities to the 6 test nodes based on theta evaluated
F6node=B(1)+testdata(pos_test,1)*B(2)+testdata(pos_test,2)*B(3); % hardcoded 2d
q = ones(length(pos_test),1)./(ones(length(pos_test),1)+exp(-F6node))


interval = 0.25;
tau = 0.5;
B1 = [B(1)-tau*.9:tau*0.3:B(1)+tau*.6];
B2 = [B(2)-2:interval:B(2)+1.5];
B3 = [B(3)-2:interval:B(3)+1.5];

%training_loss_grid;
%route_cost_grid;
%adding_the_two_costs;
%legacy end



